
// this is autogenerated file, do not edit it.
#include "ficus/ficus.h"


/*
    Linear system solution using QR factorization
    Converted from OpenCV implementation, here is the original copyright:
    ---
    This file is part of OpenCV project.
    It is subject to the license terms in the LICENSE file found in the top-level directory
    of this distribution and at http://opencv.org/license.html
    ---
*/
enum {
    FX_MATX_BUF = 1024,
};

#define FX_MATX_QR_EPS 1e-12

static int fx_decomp_qr(const double* A_, int_ astep, int_ m, int_ n, int_ k,
                        const double* b, int_ bstep, double* x, int_ xstep,
                        double* hFactors, double scale, double eps,
                        double* det_)
{
    int status = 1;
    double localbuf[FX_MATX_BUF], *buf = localbuf;
    double *A, *vl;
    int_ bufsize = m*n + m + n, nm = m < n ? m : n;

    if (det_) *det_ = 0.;
    if (bufsize > FX_MATX_BUF) {
        buf = (double*)fx_malloc(bufsize);
        if (!buf) return FX_SET_EXN_FAST(FX_EXN_OutOfMemError);
    }
    A = buf;
    vl = buf + m*n;
    for(int i = 0; i < m; i++)
        memcpy(A + i*n, A_ + i*astep, n*sizeof(A[0]));
    astep = n;

    if (!hFactors)
        hFactors = vl + m;

    for (int_ l = 0; l < n; l++) {
        //generate vl
        int_ vlSize = m - l;
        double vlNorm = 0.;
        for (int_ i = 0; i < vlSize; i++) {
            vl[i] = A[(l + i)*astep + l];
            vlNorm += vl[i] * vl[i];
        }
        double tmpV = vl[0];
        vl[0] = vl[0] + (vl[0] < 0. ? -1 : vl[0] > 0.)*sqrt(vlNorm);
        vlNorm = sqrt(vlNorm + vl[0] * vl[0] - tmpV*tmpV);
        for (int_ i = 0; i < vlSize; i++)
            vl[i] /= vlNorm;
        //multiply A_l*vl
        for (int_ j = l; j < n; j++) {
            double v_lA = 0.;
            for (int_ i = l; i < m; i++)
                v_lA += vl[i - l] * A[i*astep + j];

            for (int_ i = l; i < m; i++)
                A[i*astep + j] -= 2 * vl[i - l] * v_lA;
        }

        //save vl and factors
        hFactors[l] = vl[0] * vl[0];
        for (int_ i = 1; i < vlSize; i++)
            A[(l + i)*astep + l] = vl[i] / vl[0];
    }

    if (x) {
        for(int_ i = 0; i < m; i++) {
            if(b)
                memcpy(x + i*xstep, b + i*bstep, k*sizeof(x[0]));
            else {
                memset(x + i*xstep, 0, k*sizeof(x[0]));
                x[i*xstep + i] = scale;
            }
        }
        //generate new rhs
        for (int_ l = 0; l < n; l++) {
            //unpack vl
            vl[0] = (double)1;
            for (int_ j = 1; j < m - l; j++)
                vl[j] = A[(j + l)*astep + l];

            //h_l*x
            for (int_ j = 0; j < k; j++) {
                double v_lB = (double)0;
                for (int_ i = l; i < m; i++)
                  v_lB += vl[i - l] * x[i*xstep + j];

                for (int_ i = l; i < m; i++)
                    x[i*xstep + j] -= 2 * vl[i - l] * v_lB * hFactors[l];
            }
        }
        //do back substitution
        for (int_ i = n - 1; i >= 0; i--) {
            for (int_ j = n - 1; j > i; j--) {
                double f = A[i*astep + j];
                for (int_ p = 0; p < k; p++)
                    x[i*xstep + p] -= x[j*xstep + p] * f;
            }
            if (fabs(A[i*astep + i]) < eps) {
                status = 0;
                break;
            } else {
                double f = A[i*astep + i];
                for (int_ p = 0; p < k; p++)
                    x[i*xstep + p] /= f;
            }
        }
    }

    if(det_) {
        double p = nm & 1 ? -1 : 1;
        for(int_ i = 0; i < nm; i++)
            p *= A[astep*i + i];
        *det_ = p;
    }

    if(buf != localbuf)
        fx_free(buf);

    return status;
}

FX_EXTERN_C int fx_init_Array(void)
{
   int fx_status = 0;
   return fx_status;
}

FX_EXTERN_C void fx_deinit_Array(void)
{

}

